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We outline a Kohn-Sham-Dirac density-functional-theory (DFT) scheme for graphene sheets that 
treats slowly-varying inhomogeneous external potentials and electron-electron interactions on an 
equal footing. The theory is able to account for the the unusual property that the exchange- 
correlation contribution to chemical potential increases with carrier density in graphene. Conse- 
quences of this property, and advantages and disadvantages of using the DFT approach to de- 
scribe it, are discussed. The approach is illustrated by solving the Kohn-Sham-Dirac equations 
self-consistently for a model random potential describing charged point-like impurities located close 
to the graphene plane. The influence of electron-electron interactions on these non-linear screening 
calculations is discussed at length, in the light of recent experiments reporting evidence for the 
presence of electron-hole puddles in nearly-neutral graphene sheets. 
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I. INTRODUCTION 

Graphene is a newly realized two-dimensional (2D) 
electron system^ which has engendered a great deal 
of interest because of the new physics which it exhibits 
and because of its potential as a new material for elec- 
tronic technology. The agent responsible for many of 
the interesting electronic properties of graphene sheets is 
the non-Bravais honeycomb-lattice arrangement of car- 
bon atoms, which leads to a gapless semiconductor with 
valence and conduction 7r-bands. States near the Fermi 
energy of a graphene sheet are described by a massless 
Dirac equation which has chiral band states in which the 
honeycomb-sublattice pseudospin is aligned either par- 
allel to or opposite to the envelope function momen- 
tum. The Dirac-like wave equation leads to both un- 
usual electron-electron interaction effects and to unusual 
response to external potentials. 

Many new ideas that are now being explored in 
graphene electronics are still based on idealized mod- 
els which neglect disorder and electron-electron inter- 
actions, and as a consequence many of these may ulti- 
mately require qualitative and quantitative revision as 
our understanding of this material improves. In this pa- 
per we outline one approach, a Kohn-Sham-Dirac DFT 
scheme, which can be used for more realistic modelling 
of graphene sheets, including both disorder and electron- 
electron interactions. 

Because of band chirality, the role of electron-electron 
interactions in graphene sheets differs in some essential 
way d 3 * 4 * 5 -! from the role which it plays in an ordinary 2D 
electron gas. One important difference is that the con- 
tribution of exchange and correlation to the chemical po- 
tential is an increasing rather than a decreasing func- 
tion of carrier- density. As we discuss later, this property 
implies that exchange and correlation increases the ef- 
fectiveness of screening, in contrast to the usual case in 
which exchange and correlation weakens screening. This 



unusual property follows from the difference in sublat- 
tice pseudospin chirality between the Dirac model's neg- 
ative energy valence band states and its conduction band 



states , and in a uniform graphene system is readily 
accounted for by many-body perturbation theory. The 
principle merit of the DFT theory we describe is that 
it allows this physics to be accounted for in graphene 
sheets in which the carrier density is non-uniform either 
by design, as in p-n junction systems^, or as a result of 
unintended disorder sources. 

A related and complementary DFT method has re- 
cently been used by Rossi and Das SarmsP to study the 
ground-state density profile of massless Dirac fermions in 
the presence of randomly-distributed charged impurities. 
Their method differs from ours in two main respects: the 
authors of Ref. 7 have (i) approximated the kinetic energy 
functional of non-interacting massless Dirac fermions by 
means of a local-density approximation (LDA) whereas in 
the present work the kinetic energy functional is treated 
exactly via the Kohn-Sham mapping (see Sect. |TT| below); 
and (ii) neglected correlation effects, which, as it will 
be clear in Sect. II B[ partly compensate the enhanced 
screening due to exchange and Dirac-equation chirality. 
Inhomogeneous graphene systems have also been studied 
using the Thomas-Fermi approximation (LDA for the ki- 
netic energy only) by Fogler and collaborators^ . 

Our paper is organized as follows. In Section [II] we 
outline the version of DFT which is appropriate for non- 
uniform carrier density graphene sheets with static exter- 
nal potentials that are smooth enough to permit neglect 
of inter-valley scattering. Many-body effects enter this 
theory via an LDA exchange-correlation potential with a 
density-dependence precisely opposite to the one famil- 
iar from ordinary LDA-DFT theory applied to parabolic- 
band inhomogeneous electron liquids. In Section |III| we 
outline the procedure we have used to solve the theory's 
Dirac-like Kohn-Sham equations. In Section [TV| we dis- 
cuss results obtained by solving the Kohn-Sham equa- 
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tions self-consistently for an illustrative random poten- 
tial model, highlighting some strengths and weaknesses 
of this approach to many-body physics in inhomogeneous 
graphene sheets. In Section [V] we briefly mention other 
problems to which the theory outlined in this paper could 
be successfully applied, comment on the relationship be- 
tween our DFT approach and ab initio DFT applied to 
graphene, and summarize our main conclusions. 



II. MASSLESS DIRAC MODEL DFT 

We consider a system of 2D massless Dirac fermions 
which are subjected to a time-independent scalar exter- 
nal potential V ex t(f)- This model applies to graphene 
sheets when the external potential varies slowly on the 
lattice-constant length scale. In this limit the exter- 
nal potential will couple identically to the two sub- 
lattices, and hence be a pseudospin scalar, and have 
negligible inter- valley scattering, justifying an envelope 
function approach 9 which promotes the perfect crystal 
Dirac bands to envelope function Dirac operators. To 
account for electron-electron interactions in graphene 
sheets, the ultrarelativistic massless-Dirac particles must 
interact via instantaneous non-relativistic Coulomb in- 
teractions. The juxtaposition of an ultrarelativistic free- 
fermion term and a non-relativistic interaction term in 
the Hamiltonian of a graphene sheet leads to a new-type 
of ma ny-body problem. 

DFT^SEinil i s a practical approach to many-body 
physics which recognizes the impossibility of achieving 
exact results and seeks practical solutions with adequate 
accuracy. Following a familiar line of argument ! 10 * 11 * 12 * 
which we do not reproduce here, many-body exchange- 
correlation effects can be taken into account in the 
graphene many-body problem with the same formal justi- 
fications and the same types of approximation schemes as 
in standard non-relativistic DFT 10 * 11 ! 12 !. The end result 
in the case of present interest is that ground state charge 
densities and energies are determined by solving a time- 
independent Kohn-Sham-Dirac equation for a sublattice- 



pseudospin spinor $\(r) = ((/?^(r), (p^ J (r))' L ', 
[vcr • p + l a V KS (r)] $ x (r) = e x <$>\(r) . 



(1) 



Here v ~ 10 6 m/s is the bare Fermi velocity, p = —ihV ri 
cr is a 2D vector constructed with the 2x2 Pauli ma- 
trices a i and a 2 acting in pseudospin space, I a is the 
2x2 identity matrix in pseudospin space, and VKs(r) = 
Vextir) + A7hW + KcW is the Kohn-Sham (KS) po- 
tential, which is a functional of the ground-state density 
n(r). The ground-state density is obtained as a sum over 
occupied Kohn-Sham-Dirac spinors <&\(r): 



i(r) = 4 £ $l(r)$ A (r) 



A(occ) 



^4^[|^(r)| 2 + |^(r)| 2 ], (2) 

A(occ) 



where the factor 4 is due to valley and spin degenera- 
cies and {ip^\r),<j = A,B} are the pseudospin (sublat- 
tice) components of the Kohn-Sham-Dirac spinor §\(r). 
Equation ^ is a self-consistent closure relationship for 
the Kohn-Sham-Dirac equations ([!]), since the effective 
potential in Eq. (JlJ is a functional of the ground-state 
density n(r). More explicit details on the construction 
of n(r) are given below. This formalism is readily gener- 
alized to account for spin-polarizatiorP^l, or valley polar- 
ization 14 , or both. A generalization of the present the- 
ory to situations in which graphene is subjected to an 
inhomogenous magnetic field (as in magnetically-defined 
graphene quantum dots 15 ) can also be envisioned along 
the lines of e.g. Ref. [T6l 

The KS potential Vks(t) m Eq.( IJ is the sum of ex- 
ternal, Hartree, and exchange-correlation contributions. 
The Hartree potential 



A^ H (r) = J d 2 r' 



e\r ■ 



Sn(r f ) 



(3) 



where e is the average background dielectric constant 
(e = 2.5, for example, for graphene placed on Si02 with 
the other side being exposed to air) and the quantity 
Sn(r) = n(r) — no is the density measured relative to that 
of a uniform neutral graphene sheet as specified more pre- 
cisely below [see Eq. ([34])] . 

The third term in Vks(v), Kc(^) 5 is the exchange- 
correlation potential, which is formally a functional of 
the ground-state density, but known only approximately. 
In this work we employ the local-density approximation, 



V xc (r) 



,,hom 



n^n c (r) 



(4) 



where ^° m (n) is the reference exchange-correlation 
potential of a uniform 2D liquid of massless Dirac 



fermions^ with carrier density 



,,hom 



(n) is related 



to the ground-state energy per excess carrier fe xc (n) by 



.horn 



(n) 



<9[nfe xc (n) 
dn 



(5) 



The carrier density n c (r) is the density relative to that of 
a uniform neutral graph ene sheet and will be defined more 
precisely in Sect. IV A| The expression used for fe xc (n) 
depends on the zero-of-energy, which is normall}0H' cho- 
sen so that v^ m (n = 0) = 0. 

To apply the LDA-DFT formalism to graphene it is 
necessary to have convenient expressions for the excess 
exchange-correlation energy 5s xc (n), which will be pro- 
vided below in Sects. II A and |IIB| This quantity has 
been calculated at the Random Phase Approximation 
(RPA) level in Ref. [3j 



A. Exchange potential 

Because the Coulomb energy and the Dirac band en- 
ergy scale in the same way with length, we can write the 
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first-order exchange contribution to Se xc (n) as 

fe x (n) = e F a gT F(A) . (6) 

Here e F = sgn(n) hvk F is the Fermi energy, where k F = 
(47r|n|/g) 1//2 is the Fermi wave vector corresponding to 
an electron (hole) density n above (below) the neutrality 
point and g = g s g v = 4 accounts for spin and valley 
degeneracy. The quantity a gr is defined by a gr = g x 
e 2 / (eHv) = gx a ee , where a ee is graphene's fine structure 
constant. The ultraviolet cut-off A in Eq. Q is defined 
by A = k max/^F? where /c ma x should be assigned a value 
corresponding to the wavevector range over which the 
continuum model describes graphene. For definiteness 



we take & max to be such that 



irk n 



A 



(7) 



where Ao = ?>V?>Oq/2 ~ 0.052 nm 2 is the area of the 
unit cell in the honeycomb lattice, (ao — 1.42 A is the 
Carbon-Carbon distance) and n is a dimensionless num- 
ber n G (0, 1]. The optimal value of n would have to be 
determined by a lattice-model correlation energy calcu- 
lation. From another point of view 77, the Dirac velocity 
and the dielectric constant e are coupled parameters 
of the Dirac model for graphene which should be fixed by 
comparison of the model's predictions with experiment. 
For typical graphene-system densities, the dependence of 
the exchange-correlation potential on n is weak enough 
that we can arbitrarily choose n = 1 with some confi- 
dence. Given a value of 77, the dependence of A on density 
is given by 



AH = Vgv 



V\ n \Ao ' 



(8) 



The exchange potential corresponding to Eq. ([6| is 
given by 



dn 
dF 

dA 



2 
OA 

dn 



where 



n 



OA 

dn 



-A . 



(9) 



(10) 



We have chosen the following simple formula for F(A) to 
parametrize the data in Ref. [3j 



(11) 



F(A) = — ln(A) + 7^— 

y J 6g K } ^ 1 + b e A c 



where the first term, which is the leading contribution 
in the limit A ^> 1, has been calculated analytically in 
Ref. [3j This term is largely responsible for the quasipar- 
ticle velocity enhancement in doped graphene sheets^'. 
The numerical constants a e , 6 e , and c e are given by 



a e = 0.0173671 
b e = 3.6642 x 10- 
c e = 1.6784 



(12) 



Eq. (11) implies that 



dF 
OA 



A c 



1 1 



(1 + b e A c -) 2 A 6g A 



(13) 



Note that for n — > the exchange potential goes to zero 
like 



v* om (n 0) oc -sgn(n)a gr y/\n\ In \n\ , (14) 
i.e. with an infinite slope. 

B. RPA correlation potential 

The RPA correlation energy data of Ref. |3|can be con- 
veniently parametrized by the following formula 



Sef PA (n) 



€ F 



+ 



--^eKr)ln(A) 

Q^g r ftc(^gr) 

1 + b c {a sr )A c "M 



(15) 



where 



' a c (a gr ) = -1/(63.0963 + 57.351226 a gr ) 
b c (a g r) = (7.75095 - 0.08371 g^ 61167 ) x lO" 7 
c c (a gr ) = 1.527 + 0.0239 a gr - 0.001201 a 2 



(16) 



and 



£(a gr ) 



^ r+oc 

= 2 Jo 



dx 



(1 + x 2 ) 2 {\/l + x 2 + 7ra gr /8) 



.(17) 



Once again, the logarithmic contribution in Eq. ( 15 ) rep- 



resents the leading term in the limit A ^> 1 and was 
calculated analytically in Ref. [3j 

Note that we can write Eq. (Il5|) in the form: 



fe* PA = e F a 2 ST G agt (A) 



with 



G agr (A) = -^ln(A) 



a c (oigr) 



1 + 6 e (a gr ) Antes') 



(18) 



(19) 



Following the same procedure highlighted above for the 
exchange contribution, one easily finds the correlation 
contribution to ^xc m ( n )- The on ly necessary input to 
calculate this contribution is 



0G o 



A c = t j(a gI ) 1 
6g A 



dA (l + b c A c c)2 A 

In the limit rnOwe find 

v^ om (n -> 0) oc sgn(n)a 2 r t(a gI )^y\n\ In \n\ 



(20) 



(21) 



A plot of the exchange-correlation potential as a function 
of the density n is given in Fig. [I] For the sake of compar- 
ison, in Fig. [l] we also have plotted the quantum Monte 
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Carlo exchange-correlation potential of the parabolic- 
band 2D electron gas^, after having antisymmetrized 
it for n < 0. We can clearly see from this plot how 
the density-dependence of the exchange-correlation po- 
tential of a uniform 2D liquid of massless Dirac fermions 
is precisely opposite to the one familiar from the ordinary 
LDA for parabolic-band inhomogeneous electron liquids. 
While the latter is negative for positive density, favor- 
ing inhomogeneous densities, the former increases the 
energy cost of density increases, favoring more homoge- 
neous densities and enhancing screening. It is also ap- 
parent from this figure that the density dependence of 
the exchange-correlation potential can in some circum- 
stances lead to effects which can give the appearance of 
a gap in the graphene sheet's Dirac bands. 




n (10 12 cm" 2 ) 



III. KOHN-SHAM-DIRAC EQUATION 
SOLUTIONS: PLANE- WAVE METHOD 

In this Section we discuss Kohn-Sham-Dirac equation 
solutions based on a supercell method and plane-wave 
expansions. We consider massless Dirac fermions in a 
2D (square) box of size L x L with periodic boundary 
conditions. In this case the Kohn-Sham-Dirac equations 
([!]) can be conveniently solved by expanding the spinors 
&\(r) in a plane- wave basis. We discretize real space: 
r -> rij = (xi,yj), Xi = i5x, yj = j5y with i = l...N x 
and j = l...N y . Here Sx x N x = Sy x N y = L. Fourier 
transforms f(k) of real-space functions f(r) are calcu- 
lated by means of a standard fast-Fourier-transform algo- 
rithm 18 that allows us to compute / on the set of discrete 
wave vectors fc™, 
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10 



2tt 

T 



(22) 



with -N x /2 < n Xii < N x /2 and -N y /2 < n y j < N y /2 
(or, equivalently, < n x ^ < N x and < n y j < N y ), 
that belong to the Bravais lattice of the discretized box. 
The definition of the Fourier transform that we use is the 
following: 



f(k) = J d 2 r f(r) e~ ikr 



(23) 



After discretization f(r) -> fa = /(r^-), fa = f(kij) 
with 



FIG. 1: (Color online) Top panel: the exchange and RPA cor- 
relation potentials, v^ om (n) [(black) solid line)] and i;J om (n) 
[(blue) dashed line)], (in meV units) as functions of the den- 
sity n (in units of 10 12 cm -2 ) for a ee = 0.5. Note how for 
n — > both potentials have an infinite slope. The (magenta) 
dash-dotted line represents the full exchange-correlation po- 
tential, v^ m (n) = ^ om (n) + ^ om (n). The (green) dotted line 
is the quantum Monte Carlo exchange-correlation potential of 
a standard parabolic-band 2D electron gas^. For convenience 
we have chosen parameters corresponding to a 2D electron 
gas on a background with dielectric costant 4 and with band 
mass 0.067 m, m being the electron mass in vacuum. Bottom 
panel: the full exchange potential [(black) solid line)] is com- 
pared with its ln-only approximation, [(blue) dashed line] , i. e. 
retaining only the first term in Eq. ( 11 ). 



N X ~lNy-l 



/ v / v J nrn 

n=0 m=0 



and 



N X ~l Ny-1 



(24) 



nm e ik ^-r nm _ ( 2 5) 



In all the numerical calculations reported on below we 
use L as unit of length, 2ttH/L as the unit of momentum, 
and hv/L as the unit of energy. In what follows we also 
set h = 1. 

In momentum space Eq. reads 

^(/c|[^-p + I^Ks(r)]|/c / )^A(/c / )=^^AW . (26) 
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Here A labels the eigenvalues of the Kohn-Sham-Dirac 
matrix = (k\[va • p + I a V K s(r)]\k') . The matrix 

elements of the kinetic Hamiltonian are given by 



(fe| va - p \k') = va ■ k' '5k,k' • 



(27) 



We employ a momentum space cut-off k Xy i,k y j G 
[— k c ,+k c ] which does not exceed the Brillouin-zone 
boundary defined by our real-space discretization: k c < 
7r/5x,7r/5y. k c defines the range of momenta used in the 
expansion of the Hamiltonian W^^? and thus defines its 
dimension d^: 



4 = 2x 2 x 



Lk c 



1 



(28) 



The factor of 2 here is due to the sublattice pseudospin 
degree of freedom. As already stated above, real spin 
and valley degrees of freedom enter our calculations only 
through the trivial degeneracy factors they imply. Given 
a value of k c the Kohn-Sham-Dirac matrix 7^ |P has dn 
eigenvalues, labeled by the discrete index A = 1, . . . , ^h- 



IV. NON-LINEAR SCREENING OF COULOMB 
IMPURITIES 

As an illustration we apply the LDA-DFT method 
described above to study the non-linear screening of 
A^imp > 1 point-like impurities with charge Ze (Z can 
be either positive and negative and e > in this work) 
located at random positions on a plane at a distance d 
from the 2D chiral electron gas (CEG) plane. The ap- 
proximately linear dependence of conductivity on car- 
rier density in graphene sheets suggest d 19 ! 20 ! that nearby 
charged impurities are the dominant disorder source in 
most current graphene samples. 



A. Constructing the KS potential and the 
ground-state density 

We assume that the 2D CEG has a spatially averaged 
7r-electron density 

2 



n 



(29) 



Here 2/Ao is the density of a neutral graphene sheet and 
n c is the spatially averaged carrier density, which can be 
positive or negative and controlled by gate voltageg^ |21 | 22 [ 
In what follows we write n c = 4Q/L 2 , where Q is the 
number of carriers per spin and valley in our supercell. 
Because of the role played by gate voltages in experi- 
ment, there is no reason to impose a charge-neutrality 
relationship between the number of impurities A^ mp and 

Q. 

The external potential V ext (r) is given by: 



iVir 



^ext(r) 



Ze 2 



\ e^\r-R % \ 



d 2 



(30) 



where Ri are random positions in the supercell. For sim- 
plicity, all charges have been taken to have the same Z in 
Eq. ( |30|). T he matrix elements of the disorder potential 
in Eq. (|3Q[) are given by 

(fe|^ext(r)|fe , > = V^ t (k-k')T^ v (k-k') , (31) 

where V ext (q) = — 27rZe 2 exp (— qd)/{eq) is the Fourier 
transform of the potential created by a single impurity 
and 



*^"imp(^ 



N- 

_l J "imp 

fe') = ^ E e- i(k - k "> Ri 

i=l 



(32) 



is a geometric form factor that depends only on the po- 
sitions of the impurities. The impurity charges are repli- 
cated in each supercell and the total potential V ext (r) 
therefore has the supercell periodicity. We set V ext (k = 
k') = 0, thereby choosing the zero of energy at the Dirac- 
point energy in the spatially averaged external potential. 

The ground-state density profile n(r) in the external 
potential given by Eq. (30) is computed from Eq. ^ 
by summing over A = 1, . . . , A max , where the KS en- 
ergy levels are arranged in ascending order, e\ < • • • < 
£ A max < ' • • < £d H • Since half of the system's 7r-orbitals 
are occupied in a neutral graphene sheet, A max is related 
to the average 7r-electron density of the graphene sheet 
n = 4(d H /2 + Q)/L 2 by 



dy 



Q ■ 



(33) 



Note that this implies the following relationship between 
the momentum- space cutoff k c and the area of the sys- 
tem L 2 in units of A : L 2 /A = 2 [2L/c c /(2tt) + l] 2 . 
In our self-consistent numerical calculations we evaluate 
only the deviation of the density from its average value 
in the supercell: 



Sn(r) — n(r) — no 



(34) 



The corresponding quantity in momentum space 5n(k) 
is given by 8n(k) = n(k) — no#fc,o- Note that 5n(r) is 
charge neutral, i.e. 5n(k = 0) = 0. The matrix elements 
of the Hartree term in the Kohn-Sham-Dirac equation 
are given by 

27re 2 

<fc|AVk(r)|fc'> = Sn(k - k') . (35) 

The matrix elements of the exchange-correlation po- 
tential can be calculated numerically from 

<fc|V xc (r)|fc') ^jjij^ y xc(r) e -*(*-*'>-r , (36) 

where V xc (r) is given by Eq. Q with the carrier density 

n c (r) = n(r) = Sn(r) + . (37) 
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B. Numerical results 

In this Section we report some illustrative numerical 
results that we have obtained applying the LDA-DFT 
method described above. All the numerical results pre- 
sented in this work were obtained with 77 = 1 [see Eq. ^ 
for the definition of 77]. 

In Fig. [2] we illustrate the real-space density profile 
5n(r) of a neutral-on-average (Q = 0) 2D CEG sub- 
jected to the external potential of Af imp = 40 impuri- 
ties with Z = +1 located at a distance d = 0.1 L from 
the graphene plane [the corresponding external potential 
V e xt(r) is illustrated in the top left panel of Fig. [2]. In this 
particular simulation we have used a ee = 0.5 and k c = 
(2tt/L) x 10, which corresponds to an effective square size 
L 2 = 882 Ao ~ 46 nm 2 . The charges are therefore sepa- 
rated from the graphene layer by d ~ 0.7 nm. This model 
is motivated by growing experimental evidence that the 
dominant source of disorder in most graphene samples 
is external charges, probably located in the nearby sub- 
strate. 

In Fig. [2] we have reported: (i) the non-interacting 
Dirac electron density profile, which is obtained by 
setting the Hartree and exchange-correlation potentials 
in the Kohn-Sham-Dirac Hamiltonian to zero; (ii) the 
"Hartree-only" density profile, which is obtained by 
solving the Kohn-Sham-Dirac equations self-consistently 
with V xc (r) = 0; and (iii) the "full" density profile, 
which includes both Hartree and exchange-correlation ef- 
fects. The self-consistent calculations are iterated until 
the Kohn-Sham potential is converged to a relative pre- 
cision of ~ 10 -3 . 

Electron- hole puddles, similar to those observed in 
Refs. I21|22( are evident in all these plots, although 
there are qualitative and quantitative differences be- 
tween the non-interacting density profile and the ones 
that include electron-electron interactions (the experi- 
mental observation that the spatial pattern of electron- 
hole bubbles is not correlated with the topography o f the 
graphene sheets^, is consistent with the inferenc e 1 19 * 20 * 
from conductivity- us. -carrier density data that remote 
charges rather than sheet corrugations dominate disor- 
der). To begin with, we note how the inclusion of the 
Hartree term has the (expected) effect of reducing the 
amplitude of the spatial fluctuations of 5n(r) quite dra- 
matically, by approximately a factor of two in these non- 
linear screening calculations. It is interesting to compare 
this reduction factor with what would be expected in a 
linear screening approximation. Neutral graphene has 
the unusual property that its static dielectric function 
e(q) neither diverges as wavevector q goes to zero, as it 
would in a 2D metal, nor approaches 1, as it would in a 
2D semiconductor. Instead 

27re 2 

e(q) = 1 ~ — Xpp(q) (38) 

approaches a constant because the polarization function 
Xpp(q) ( or proper density-density response function^) 



has a non-analytic linear dependence on q due to inter- 
band transitions with vanishing energy denominators. 
In the Hartree approximation [Xpp(q) —> X^°\q)i where 
X^ \q) is the non-interacting polarization functiorffl see 
Sect. 5.3.1 of Ref . H2l for more details] 

e(q) 1 + I gctee ~ 1-78 (39) 

for the value of a ee used in our calculations. When ex- 
change and correlations corrections are included in e(q) 
increases by a small fraction, enhancing screeening. The 
influence of interactions on the non-linear screening cal- 
culations summarized in Fig. [2] is therefore (perhaps sur- 
prisingly) broadly consistent with expectations based on 
linear screening theory - even at a semi-quantitative level. 
Qualitative non-linear effects do however appear in some 
details, as we now explain. 

In Fig. [3] we examine the induced carrier density in 
more detail by plotting 5n(r) as a function of x for a fixed 
value of y. Here we see clearly that V xc (r) tends to cause 
the density to vary less rapidly in those spatial regions 
at which the carrier-density changes sign. The origin of 
this behavior in our calculations is that the exchange- 
correlation potential increases especially rapidly with 
density in these regions. This aspect of the induced den- 
sity profile is similar to the behavior which would be pro- 
duced by an energy-gap of ~ 0.1 eV in the graphene 
bands (see Fig. [TJ. The rapid change in exchange- 
correlation potential with density alters the statistical 
distribution of density- values in a disordered sample, as 
studied in some detail using a Thomas-Fermi approxi- 
mation for the non-interacting kinetic energy functional 
(and including local-density-approximation exchange) by 
Rossi and Das Sarma in a recent paper^. Thomas-Fermi 
theory is formally based on a gradient expansion of the 
total energy density (see e.g. Sect. 7.3.1 in Ref. H2j) . 
When applied to graphene, assuming that the typical 
length scale for density variations in the 2D CES is the 
inverse of the Thomas- Fermi screening wavevector /ctf = 
27re 2 h , (e-F) /e [here v{ey) = gsF / '(^ttv 2 ) is the density-of- 
states at the Fermi energy] , Thomas- Fermi theory can be 
viewed as an expansion in powers of &tf/^f = goL ee . As 
emphasized by Fogler and collaborators 8 , this parameter 
is not small when the value used for a ee is in the range 
~ 0.5 thought to be appropriate for graphene on Si02- 
In our approach we avoid a local-density-approximation 
for the non-interacting kinetic energy functional by solv- 
ing microscopic Kohn-Sham-Dirac equations (this is the 
idea behind the Kohn-Sham mapping 10 ). We cannot 
avoid the local-density approximation for the exchange- 
correlation potential however [Eq. Q], and it must be 
acknowledged that this is a defect of our theory, and one 
that is not easily remedied. The situation is similar to 
that in standard DFT applications, in which the local 
density- approximation is not rigorously valid on atomic 
length scales. It has nevertheless been possible to rem- 
edy defects of the local-density-approximation in many 
circumstances by using modified functionals, for exam- 
ple generalized-gradient approximations, which are often 
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FIG. 2: (Color online) Top left panel: a color plot of the ex- 
ternal potential V e ^t(r) (in units of hv/L) as a function of 
x/L and y/L. The system parameters are N x = N y = 128, 
fee = (2tt/L) x 10, iVi m p = 40, Z = +1, a ee = 0.5, Q = 0, 
and d/L = 0.1. The small circles represent the positions of 
the impurities for a particular realization of disorder. Top 
right panel: a color plot of the corresponding non-interacting 
ground-state density profile 5n(r) (in units of 1/L 2 ) as a 
function of x/L and y/L. Bottom left panel: Hartree-only 
ground-state density profile. Bottom right panel: same as in 
the bottom left panel but with the addition of the exchange 
and RPA correlation potential. 



semi-phenomenological in character. Our expectation is 
that the LDA for exchange and correlation in graphene 
will improve accuracy compared to Thomas-Fermi ap- 
proximation theories in which the band energy is also 
approximated using an LDA. In addition, it will likely 
prove possible to compensate for the main-defects of the 
exchange-correlation LDA by using modified exchange- 
correlation energy functionals which are informed by 
comparisons between theory and experiment. 

In Figs. [4]and[5]we report results similar to those pre- 
sented in Figs. [2]and[3j but for a separate realization of 
the random charged impurity distribution and a smaller 
separation between the impurity plane and the graphene 
plane, d = 0.05 L. When the impurities are closer to 
the graphene plane the role of the exchange-correlation 
potential seems to become less important. Conversely, 
for larger d exchange and correlation effects increase in 
importance. Because of the peculiar response of Dirac 
fermions, quite localized charge distributions can be in- 
duced by disorder potential features, even when those 
features are weak. Indeed we find that for large sepa- 
rations between the graphene and impurity planes, the 
Kohn-Sham-Dirac equations do not always converge, in- 
dicating the possible importance in some circumstances 
of correlation effects which cannot be captured by the KS 
LDA theory. 

Finally in Figs. [6][7| we illustrate the dependence of 
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FIG. 3: (Color online). A one-dimensional plot of 5n(r) cor- 
responding to the data in Fig. [5] as a function of x/L for 
y/L — 0.5. The circles label the non-interacting result, the 
squares label the Hartree-only self- consistent result, and the 
triangles label the full self-consistent result. 




FIG. 4: (Color online) Same as in Fig. [2] but for a different 
distribution of charges and for d/L = 0.05 instead of d/L = 
0.1. 



Sn(r) on Q, i.e. on a gate potentials which move the 
average density- away from the Dirac-point. Because 
of the unavoidable presence of external charges in any 
graphene sheet environment, this is actually the generic 
case. Special neutral sheet properties, like those re- 
ferred to below in the single impurity case, will be dif- 
ficult to realize experimentally. Fig. [6] shows the ex- 
ternal potential created by a particular distribution of 
^imp = 40 random charges, different again from the 
distributions used in producing Figs. [2p) and the cor- 
responding ground-state density profile 5n(r) calculated 
for Q = 0. The data in Fig. [6] refer to a system with 
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FIG. 5: (Color online). A one-dimensional plot of 5n(r) cor- 
responding to the data in Fig. [4] as a function of x/L for 
y/L = 0.5. The color coding is the same as in Fig. [3] 
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FIG. 6: (Color online) Same as in Figs. [2] and HI but for a 
different distribution of impurities, k c — [2n/L)x 15, and 
d/L = 0.07. 
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FIG. 7: (Color online) Illustrating the Q-dependence of Sn(r). 
Top panel: A one-dimensional plot of the non-interacting den- 
sity profile Sn(r) corresponding to the external potential in 
the top left panel of Fig. 6] as a function of x/L for y/L = 0.5. 
The (blue) circles label the result for Q = 0, the (red) squares 
label the result for Q = 10, the (green) diamonds label the 
result for Q = 20, the (cyan) triangles up label the result for 
Q = 30, and the (yellow) triangles down label the result for 
Q — 40. Bottom panel: same as in the top panel but for the 
full self-consistent density profile. 



square size L 2 = 1922 A$ ~ 100 nm 2 . We then calcu- 
late Sn(r) for the same distribution of impurities but for 
Q = 10, 20, 30 and 40. The results of these simulations 
are shown and compared in Fig. [7| From this figure we 
clearly see that increasing the average density of the sys- 
tem, increases the amplitude of the density fluctuations 
substantially when electron-electron interactions are ne- 
glected (see top panel in Fig. [7]) . When electron-electron 
interactions are included (see bottom panel in Fig. [7]), 
this effect still occurs but Sn(r) seems to saturate with 
increasing Q. Of course, the carrier density fluctuation 
decreases in a relative sense with increasing Q. 

We conclude this section, by reporting results for the 



single-impurity case. The calculation of the density dis- 
tribution of 2D non- interacting massless Dirac fermions 
in the presence of a single Coulombic impurity placed at 
the origin Ri = of the graphene pla ne (d = 0) has re- 
cently received a great deal of attentio^ 24 ' 25 | 26 | 27 | 28 | 29 | 3Q | . 
The analytical analyzes reported in these works shows 
the existence of (at least) two different regimes: (i) a 
regime termed "subcritical" , for Za ee < 1/2, in which 
the screening density Sn(r) is localized on the impu- 
rity, Sn(r) ex S(r) and (ii) a regime termed "supercriti- 
cal", for Za ee > 1/2, in which the screening density ex- 
hibits a power-law tail Sn(r) ~ 1/r 2 at large distances. 
It is important to understand how these results are al- 
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tered by the electron-electron interactions present in real 
graphene planes. The situation in graphene sheets is 
in this sense very different from standard semiconductor 
shallow-impurity problems, especially so when the Fermi 
level lies at the Dirac point. In the standard problem, 
it is a good approximation to truncate the Hamiltonian 
to a single band. Interactions then play no role since, a 
single-hole or single-electron trapped by a charged impu- 
rity does not interact with itself. In graphene, on the 
ohter hand, both conduction and valence bands must 
be retained and the single-impurity problem is really a 
many-body problem. 

The method used here to solve the Kohn-Sham-Dirac 
equations, in which we project onto a plane- wave basis, is 
not optimized for the study the single impurity problem 
because it does not take advantage of its circular sym- 
metry. Nonetheless, in Fig. [5] we present some numerical 
results for the density distribution of a 2D CEG in the 
presence of a single impurity placed at the center of the 
sample (x = L/2,y = L/2) and right on the graphene 
plane. In particular, we show a ID plot of Sn(r) as a func- 
tion of x/L for y/L = 0.5. These density profiles corre- 
spond to a Z = +1 impurity in a Dirac sea with a ee = 0.5 
and Q = 0. In the two simulation results reported in this 
figure we have used k c = (2tt/L) x 15, which corresponds 
to an effective square size of L 2 = 1922 Ao ~ 100 nm 2 
and k c = (2tt/L) x 20, which corresponds to an effective 
square size of L 2 — 3362 A§ ~ 175 nm 2 . Comparing 
the results in the top [k c = (2tt/L) x 15] and bottom 
[fc c = (2tt/L) x 20] panels we can clearly see how they are 
compatible with a completely localized screening density 
with a J-function shape, the finite-width of Sn(r) being 
solely due to our momentum-space cutoff. 

Finally, in Fig. [9] we show how 5n(r) behaves quite 
differently in the two cases a ee = 0.1 and a ee = 1.0. 
Indeed, the non-interacting density seems to possess a 
long-range tail for a ee = 1.0. When electron-electron in- 
teractions are taken into account though, it seems that 
the behavior of Sn(r) is quite similar in both cases. This 
is in agreement with the findings of Ref. 30, in which 
the authors have shown that when electron-electron in- 
teractions are taken into account at the Hartree level, 
a Z = +1 impurity always remains in the subcritical 
regime. 



V. SUMMARY AND DISCUSSION 

When inter-valley scattering is weak, doped and gated 
graphene sheets can be described using an envelope- 
function Hamiltonian with a new sublattice pseudospin 
degree-of-freedom, an ultrarelativistic massless-Dirac 
free-fermion term, a pseudospin scalar disorder poten- 
tial, and a non-relativistic instantaneous Coulombic in- 
teraction term. There is considerable evidence from ex- 
periment that this simplified description of a honeycomb 
lattice of Carbon atoms is usually a valid starting point 
for theories of those observables that depend solely on 
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FIG. 8: (Color online). One-dimensional plots of 5n(r) as a 
function of x/L for y/L = 0.5 for a single impurity with Z = 
+ 1 located at x = y = L/2. Here d = 0.0 and a ee = 0.5. Top 
panel: numerical results for k c = (2tv / L) x 15. Bottom panel: 
numerical results for k c = (2tt/L) x20. The (blue) circles label 
the non-interacting result, the (green) squares label the self- 
consistent Hartree-only result, and the (red) triangles label 
the full self-consistent result. 



the electronic properties of 7r-electrons near the graphene 
Dirac point. Although the use of this model simplifies 
the physics considerably it still leaves us with a many- 
body problem without translational invariance which we 
do not know how to solve. 

A common strategy in piecing together the physics of 
disordered interacting-fermion problems is to solve mod- 
els in which interactions are neglected, appealing perhaps 
to Fermi-liquid-theory concepts, and to solve problems in 
which disorder is neglected, hoping that it is sufficiently 
weak to be unimportant for some observations. We antic- 
ipate that this divide and conquer approach will often fail 
in graphene. With this motivation, we have presented in 
this paper a Dirac-Kohn-Sham density-functional-theory 
scheme for graphene sheets, which treats interactions and 
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FIG. 9: (Color online). One-dimensional plots of Sn(r) as a 
function of x/L for y/L — 0.5 for one impurity with Z = +1 
located dX x = y = L/2. Here d = 0.0 and /c c = (2ir/L) x 
15. Top panel: non-interacting results. Bottom panel: full 
results. In each panel the (red) triangles label the results for 
a ee = 0.1, the (blue) squares label the results for a ee = 0.5, 
while the (green) diamonds label the results for a ee = 1.0. 



smooth inhomogeneous external potentials on an equal 
footing. Although it is formally an exact solution of the 
graphene many-body problem, it relies in practice on ap- 
proximate exchange-correlation functionals. 

The best approximation available for the graphene 
problem at present is the local-density-approximation 
(LDA) for the exchange-correlation potential. In this 
paper we have provided convenient parametrizations of 
the exchange and correlation energies of uniform-density 
graphene systems based on random-phase-approximation 
many-body calculations. These results can be used to 
take account not only of density variations in a disor- 
dered graphene sheet, but also of changes in the sheets 
dielectric environment which alters the coupling constant 
which appears in the Dirac model for graphene. 

We believe that the exchange and correlation ef- 



fects captured by our DFT theory will be important 
for many-qualitative aspects of grahene electronic struc- 
ture. In graphene the dependence of the LDA exchange- 
correlation potential on density is opposite to that of nor- 
mal 2D or 3D electron systems. As explained in detail in 
Ref. [U the origin of this behavior is in the interplay be- 
tween Dirac-model free-fermion pseudospin-chirality and 
Coulomb interactions; when the carrier density is zero in 
a graphene sheet the pseudospin-chirality polarization is 
maximized and this leads to lower interaction energies. 

It is important to contrast the DFT scheme outlined in 
this paper with normal microscopic DFT applied to the 
carbon atoms of a graphene sheet. The fully microscopic 
DFT deals with all the Carbon atom orbitals, includ- 
ing the sp 2 bonding and anti-bonding orbitals, which are 
away from the Fermi level and neglected in the Dirac- 
model, and can be used for example^ to calculate the 
electron-phonon coupling in a graphene sheet from first 
principles. Microscopic DFT also provides an ab ini- 
tio estimate of the massless Dirac velocity, which is a 
phenomenological parameter of the Dirac-model theory. 
The advantages of using the present DFT scheme for 
some 7r-orbital properties of graphene sheets are made 
clear by observing that microscopic DFT, in which the 
exchange-correlation potential is based on the proper- 
ties of a uniform 3D electron gas, fails to capture the 
anomalous sign of the density- derivative of graphene's 
exchange-correlation potential. From a microscopic point 
of view this anomalous sign is a combined consequence of 
the peculiarities of Dirac bands and non-local exchange 
and correlation effects captured by the uniform-density 
Dirac-model. 

In this paper we have illustrated the properties of this 
DFT description of disordered graphene sheets by con- 
centrating on the non-uniform carrier density. Although 
the Kohn-Sham orbitals which appear in this and other 
DFT scheme are formally justified only for the role they 
play in density and ground-state-energy calculations (due 
to the Hohenberg-Kohn theorem^), their physical signif- 
icance is often interpreted more liberally by associating 
the Kohn-Sham eigenvalues with quasiparticle energies. 
This pragmatic approach can fail spectacularly, as it fa- 
mously does in the estimation of common semiconduc- 
tor band gaps, but is more often quite useful in inter- 
preting spectral properties of materials. In the case of 
7r-orbital properties of disordered graphene sheets, STM 
local-density-of-states, ARPES, and optical conductivity 
spectra require interpretation. In our view it will be use- 
ful to apply the present approach as one element of an 
effort to improve understanding of what these probes tell 
us about particular graphene sheets. The fact that the 
self-energy of uniform-density graphene sheets has a large 
dependence on wavevector relative to the Dirac-pointPl, in 
addition to its dependence on wavevector and energy rel- 
ative to the Fermi surface, may help justify taking this 
liberty with the DFT formalism. 
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